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Abstract. The formation, structure and evolution of protoplanetary 
discs is considered. The formation of giant planets within the environment 
of these models is also discussed. 



1. Introduction 

The hypothesis that the planets in the solar system were formed in a flattened 
differentially rotating gaseous disc was originally proposed by Laplace (1796). 
The suggestion arises naturally in view of their orbital configuration. They are in 
near circular orbits which approximately lie in the same plane. Models in which 
planet formation occurs in gaseous discs have been the subject of much theo- 
retical interest in recent times (e.g. Lin & Papaloizou 1985, 1993 and references 
therein). 

When discussing the large scale physical properties of such a disc, it is convenient 
to make the idealization that it is infinitesimally thin or completely flat. When, 
as in the present solar system, most of the mass lies in a central star of mass 
M*, the local gravitational acceleration can be taken to be g = GM*/r 2 . A 
gaseous element at radius r then rotates in circular orbit with angular velocity 
f2, which, to a good approximation, is obtained from centrifugal balance by 
equating g = r$7 2 . This gives Kepler's law 

« 2 = ^- a) 

In this approximation, in which forces other than that due to the gravity of 
the central star are neglected, material rotates conserving its specific angular 
momentum j = r 2 Q. This then changes only slowly when viscous or other per- 
turbative forces are taken account of. 



1.1. Observations of protostellar discs 



Between 25 and 75 percent of young stellar objects in the Orion nebula appear 
to have discs (McCaughrean & Stauffer 1994). Masses of ~ 10~ 2±1 M , and 
dimensions ~ 40 ± 20 AU have been estimated (Beckwith & Sargent 1996). 
The presence of discs on the scale of astronomical units has been inferred from 
the infrared excesses observed in about half of all T Tauri stars. Their colors and 
magnitudes can be fitted by those expected for young pre-main sequence stars 
with ages ~ 10 6 yr (Strom et al. 1993). The infrared emission may be produced 
by the gravitational potential energy liberated by matter flowing inwards at a 
rate M ~ 10" 8±1 M yr _1 . The non observation of discs around older T Tauri 
stars together with these values of M suggest a disc lifetime of ~ 10 7 yr. 

1.2. Formation of discs 

It is believed that protostars and protoplanetary discs derive from interstellar 
matter contained in molecular clouds. Observations (Goodman et al. 1993) 
indicate that typical star-forming dense cores in dark molecular clouds have 
specific angular momentum j > 6 x 10 20 cm 2 s _1 . When these clouds undergo 
gravitational collapse, j is initially approximately conserved because the collapse 
is dynamical. Gas in the outer parts will not be able to fall directly to the centre 
(r = 0) if j is not zero. 

To obtain an estimate of the initial size of the disc, we consider the idealized 
situation when the pre-collapse cloud is a cold rotating sphere of mass M and 
radius R, with a single axis of rotation. Matter located on the rotation axis has 
no angular momentum so it can fall directly to the centre. An estimate of the 
time required, tff, is given by the time required to free fall from rest through 
a distance R under the initial inward gravit ational acce leration at the surface, 
g = GM/R 2 . This gives t ff = y/2R/g = yj2R A / (GM). In terms of the initial 
mean density ~p of the cloud, tff = y/3/(2irGp). For a cloud core with M = 1 M 
and R = 0.1 pc, this gives tff ~ 6 x 10 5 yr. 

While matter located on the rotation axis can move directly to the centre, matter 
in the equatorial plane at the surface of the sphere will be at the outermost 
radius Rd of a disc after collapse. Conservation of specific angular momentum 
determines Rd- Assuming the total mass, M, to be concentrated at the centre, 
the angular velocity at the outer edge of the disc will be given by Kepler's law 
such that 

r.2 f GM 

n = = ~W (2) 

Thus Rd is given in terms of the conserved quantities j and M by Rd = j 2 / (GM). 
Adopting j = 6 x 10 20 cm 2 s" 1 and M = 1 M , we find that R d ~ 180 AU. The 
characteristic dimension of such a disc is about an order of magnitude larger 
than our present solar system but is similar to those of protostellar discs now 
being observed by direct imaging. 



2. Early evolution 

The formation of a protostellar disc through the collapse of a molecular cloud 
core takes 10 5 -10 6 yr. During the early stages when the disc is still embedded 
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(class 0/1 object) and has a significant mass compared to the central star, there 
may exist strong disc winds and bipolar outflows (e.g. Reipurth et al. 1997) 
with associated magnetic fields. During this stage a hydromagnetic disc wind 
may be an important means of angular momentum removal for the system (see 
Papaloizou & Lin 1995, and references therein). 

When the mass of the disc is significant compared to that of the star, there 
may be a short period (~ 10 5 yr) of non axisymmetric global gravitational 
instability with associated outward angular momentum transport (Papaloizou 
& Savonije 1991; Heemskerk et al. 1992; Laughlin & Bodenheimer 1994; Pickett 
et al. 1998) that results in additional mass growth of the central star. This 
redistribution may occur on the dynamical timescale (a few orbits) of the outer 
part of the disc and so may be quite rapid, on the order of 10 5 yr for R = 
500 AU. The parameter governing the importance of disc self-gravity is the 
Toomre parameter, Q = M*H/ (M^r), with being the disc mass contained 
within radius r and H being the disc semi-thickness. In this review we shall take 
H to be the distance between the disc mid-plane and surface. This is usually 
a factor of 2-3 larger than the vertical extent reached by a disc particle moving 
through the disc mid-plane with the local sound speed c s . Thus H ~ (2c s )/fi. 
Typically H/r ~ 0.1 (Stapelfeldt et al. 1998) such that the condition for the 
importance of self-gravity, Q ~ 1, gives ~ 0.1M*. 

The characteristic scale associated with growing density perturbations in a 
disc undergoing gravitational instability with Q ~ 1 is ~ H, and the correspond- 
ing mass scale is M^H/r) 2 ~ M*(H/r) 3 , which is ~ 1 Mj for H/r ~ 0.1 and 
M* = 1 M , with Mj being Jupiter's mass. Gravitational instability does not 
necessarily lead to fragmentation (e.g. Pickett et al. 1998), nonetheless it has 
been proposed as a mechanism for directly forming giant planets by Cameron 
(1978) and more recently by Boss (1998). 

During the period of global gravitational instability, it is reasonable to sup- 
pose that the disc mass is quickly redistributed and reduced such that global 
gravitational stability is restored (Q > 1), after which further disc evolution oc- 
curs on a longer timescale governed by viscosity with effects due to self-gravity 
being small. 

2.1. Viscous evolution 

During this phase, the disc may attain a configuration similar to that expected 
for the minimum mass primordial solar nebula, ~ 10 _2 -10 _1 Mq. Planets 
have been proposed to form out of such a disc by a process of growth through 
planetesimal accumulation followed, in the giant planet case, by gas accretion 
(Safronov 1969; Wetherill & Stewart 1989). 

During this evolutionary phase, it is reasonable to regard the disc as an 
axisymmetric configuration in which, to a first approximation, material orbits 
in circlar Keplerian orbits. However, other weaker forces due to internal pressure, 
viscosity, or magnetic fields may also operate in the disc. These can result in 
angular momentum redistribution on a long timescale. In order to flow inwards, 
material has to transport its angular momentum outwards to matter at larger 
radii. The angular momentum transport process determines the timescale on 
which mass accretion can occur and hence the evolutionary timescale of the disc. 
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Historically, the first angular momentum transport mechanism to be considered 
was through the action of viscosity (von Weizsacker 1948). This acts through 
the friction of neighbouring sections of the disc upon each other. The inner 
regions rotate faster than the outer regions and thus viscous friction tends to 
communicate angular momentum from the inner parts of the disc to the outer 
parts. In order to result in evolution on astronomically interesting timescales, it 
is necessary to suppose that an anomalously large viscosity is produced through 
the action of some sort of turbulence. The magnitude of the viscosity is usually 
parameterized through the Shakura & Sunyaev (1973) a prescription. 
In this we suppose the kinematic viscosity coefficient v = ac^/fl ~ aH 2 Q, 
where a is a dimensionless constant which must be less than unity. Currently 
the most likely mechanism for producing turbulence is through hydromagnetic 
instabilities (Balbus & Hawley 1991) which might produce a ~ 0.01, provided 
the disc has adequate ionization. 

Gammie (1996) has proposed a model in which viscosity only operates in 
the surface layers where external sources of ionization such as cosmic rays can 
penetrate. Such layered models may have an interior dead zone of material for 
r > 0.1 AU and, although they can be considered as models with a variable a, 
they will behave somewhat differently from the standard models considered here 
which have a assumed to be constant throughout. 



3. The diffusion equation for disc evolution 

The evolution of a viscous disc is controlled by angular momentum conservation. 
The conservation of specific angular momentum may be expressed in the form ( 
eg. Papaloizou & Lin 1995) 

Di Id d 

3 = _ V • F = -~(rF r ) - jr-{F z ), (3) 
Dt r or oz 

where F = (F r , F z ) is the angular momentum flux, with the vertical coordinate 
being denoted by z, and p is the mass density. 

The process of vertical averaging applied to equation (||) in the case of a Keple- 
rian disc yields 

= W»> ( 4 ) 

ar r or 

where v r is the radial velocity in the disc and X = J^ H pdz is the surface mass 
density. The vertical average for a quantity Q is defined by 



I-h pQ- dz 

f-H P dz 



(Q) = ~H ■ ( 5 ) 



We assume F vanishes on the disc boundaries ensuring that the total angular 
momentum of the disc is conserved. We comment that this implies that angular 
momentum loss through winds or gain through mass infall is neglected. The 
radial angular momentum flux arising from viscosity is given by 
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Fr = " rV *' (6) 



and thus 



(Fr) = -r 2 (,)Sp (7) 
dr 

For a Keplerian disc the angular velocity decreases outwards so that the angular 
momentum flux is directed outwards as required in order that mass accretion 
may occur. 

Using equation (|^) in equation (||) and solving for (v r ) gives 

rS \dr J dr \ dr J 



For a Keplerian disc in which j = \/GM*r, we get 

Thus the radial velocity in the disc is determined in terms of the kinematic 
viscosity. For a steady state disc in which (i')E is constant, we find 

3(1/) 



(i; r ) 



2r 



This velocity is negative, consistent with accretion onto the central object. 

To obtain a general equation governing the disc, the velocity (||) is used together 

with the vertically averaged continuity equation 

<9S 1 d . 

_ + -- (&(„•■» =0. 

The global evolution of the disc is thus found to be governed by a single diffusion 
equation for the surface density which takes the form (Lynden-Bell & Pringle 
1974) 



<9i r dr \ dr 

According to this the characteristic evolution timescale for the disc will be the 
global diffusion timescale. The diffusion coefficient is 3(is) and, adopting a char- 
acteristic radius r, the global diffusion timescale is 

r 2 



3(u) 



For viscosity coefficient u, parameterized through the a prescription so that 
v = ac^/Q ~ aH 2 Q, we obtain 

t e n = (l/3)(r/H) 2 a-\ 
5 



For a protostellar disc of size 100 AU with H/r ~ 0.1, and a central solar 
mass, we find t e = 5 x 10 5 (0.01/a) yr. This gives lifetimes comparable to those 
estimated for discs around T Tauri stars if a ~ 10~ 3 -10~ 2 . 

Equation (^) enables the evolution of the disc to be determined if (u) is 
specified as a function of E. This can be done by solving for the vertical structure 
of the disc. 



4. Vertical structure calculations 
4.1. Basic equations 

We adopt the equation of vertical hydrostatic equilibrium: 

~ = (10) 

p oz 

and the energy equation, which states that the rate of energy removal by radia- 
tion is locally balanced by the rate of energy production by viscous dissipation: 

— = -pvQ,\ 11 
oz 4 

where J- is the radiative flux of energy through a surface of constant z which is 
given by: 

-16aT 3 dT 
3kp dz 

Here P is the pressure, T is the temperature, k is the opacity, which in general 
depends on both p and T, and a is the Stefan-Boltzmann constant. 
To close the system of equations, we relate P, p and T through the equation of 
state of an ideal gas: 

P = ^, (13) 
pniH 

where k is the Boltzmann constant, p is the mean molecular weight and rriH 
is the mass of the hydrogen atom. Since the main component of protostellar 
disks at the temperatures we consider is molecular hydrogen, we take p = 2. We 
denote the isothermal sound speed by c s {c 2 s = P/p). 

As above, we adopt the a-parametrization of Shakura & Sunyaev (1973), 
so that the kinematic viscosity is written v = a(? s j£l. In general, a is a function 
of both r and z. However, we shall limit our calculations presented below to 
cases with constant a. With this formalism, equation (O) becomes: 



dT 9 

oz 4 

Boundary conditions We have to solve three first order ordinary differential 
equations for the three variables P (or equivalently p), and T as a function 
of z at a given radius r. Accordingly, we need three boundary conditions at each 
r. We shall denote with a subscript s values at the disk surface. 
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The first boundary condition is obtained by integrating equation (11) over 
z between —H and H. Since by symmetry T{z = 0) = 0, this gives: 

T s = ^-M st n 2 , (15) 

where we have defined M st = 37r(i/)£. If the disk were in a steady state, M st 
would not vary with r and would be the constant accretion rate through the 
disk. In general however, this quantity does depend on r. 

The second boundary condition is determined by the fact that very close 
to the surface of the disc, since the optical depth r a £, above the disc is small, we 
have: 

Ps = 9s T ab/ K s • 

This condition is familiar in stellar structure (e.g. Schwarzschild 1958). 
Using g s = Q 2 H, we thus obtain 

P s = ^. (16) 

Provided Tjj « 1, the results do not depend on the value of T a b we choose (see 
below). 

The third boundary condition we use is 
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9akT s tt 3 

Here the disc is assumed immersed in a medium with background temperature 
If,. The surface opacity k s in general depends on T s and p s and we have used 
c 2 s = kT/ (/zm#). 

The boundary condition ([H]) is the same as that used by Levermore <fe 
Pomraning (1981) in the Eddington approximation (their eq. [56] with 7 = 1/2). 
In the simple case when T b = 0, and the surface dissipation term involving a 
is set to zero, it simply relates the disc surface temperature to the emergent 
radiation flux. 



Disc models At a given radius r and for a given values of the parameters 
M s t and a, we solve equations (|l0[), ( |l2| ) and (|14|) with the boundary condi- 
tions (|l5|), (p"6|) and ( |T7| ) to find the dependence of the state variables on z. 
The opacity is taken from Bell & Lin (1994). This has contributions from dust 
grains, molecules, atoms and ions. It is written in the form k = Kip a T b where 
Ki, a and b vary with temperature. 

The equations are integrated using a fifth-order Runge-Kutta method with 
adaptive step length (Press et al. 1992). For a specified M s t, we determine 
the value of H, the vertical height of the disc surface. This is done iteratively. 
Starting from an estimated value of H, after satisfying the surface boundary 
conditions, the equations are integrated down to the mid-plane z = 0. The 
condition that J- = at z = will not in general be satisfied. An iteration 
procedure (e.g. the Newton-Raphson method) is used to adjust value of H until 
jT = 0atz = 0toa specified accuracy. 
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An important point to note is that as well as finding the disc structure, we 
also determine the surface density E for a given M s t = 37r(f)S. Thus a (z^),£ 
relation is derived. 

In the calculations presented here, we have taken the optical depth of the 
atmosphere above the disk surface T a b = 10 -2 and a background temperature 
Tb = 10 K. In these calculations the temperatures are lower than about 4,000 K, 
so that, at the densities we consider, hydrogen is not dissociated and the mean 
molecular weight /i = 2. 

In the optically thick regions of the disk, the value of H is independent of the 
value of T a b w e choose. However, this is not the case in optically thin regions 
where we find that, as expected, the smaller r a b, the larger H. However, this 
dependence of H on r a b has no physical significance, since in all cases, the mass 
is concentrated towards the disk mid-plane in a layer with thickness independent 
of T ah . 

For example, for M st = 10~ 8 M© yr -1 and a = 10~ 2 , we find, at r = 100 AU, 
that H/r = 0.08 and 0.24 for T a b = 10 -2 and 10 -5 , respectively. However, in 
both cases, the surface density, the optical thickness and the mid-plane temper- 
ature are the same. Only the mid-plane pressure varies slightly (by about 30%) 
between these cases. 

In Figures |l]a-c, we plot H/r, S and the mid-plane temperature T m versus 
r for M s t between 10 -9 and 10 -6 M Q /year (assuming this quantity is the same 
at all radii, i.e. the disk is in a steady state) and for illustrative purposes we 
have adopted a = 10~ 3 . 

Figure ||a indicates that the outer parts of the disk are shielded from the 
radiation of the central star by the inner parts. This is in agreement with the 
results of Lin & Papaloizou (1980) and Bell et al. (1997). For a = 10" 3 , the 
radius beyond which the disk is not illuminated by the central star varies from 
0.2 AU to about 3 AU when M st goes from 10" 9 to 10" 6 M Q yr _1 . These values 
of the radius move to 0.1 and 2 AU when a = 10~ 2 . Since reprocessing of 
the stellar radiation by the disk is not an important heating factor below these 
radii, this process will in general not be important in protostellar disks. We 
note that this result is independent of the value of T a b we have taken. Indeed, as 
we pointed out above, only the thickness of the optically thin parts of the disk 
(which do not reprocess any radiation) gets larger when r a b is decreased. 

The values of H/r, £ and T m we get are qualitatively similar to those 
obtained by Lin & Papaloizou (1980), who adopted a prescription for viscosity 
based on convection, and Bell et al. (1997). Our values of H/r are somewhat 
larger though, since H is measured from the disk mid-plane to the surface such 
that T a b is small, and not 2/3 as usually assumed. However, as we commented 
above, this has no effect on the other physical quantities. We also recall that H, 
as defined here, is about 2-3 times larger than c s /£l, with c s being the mid-plane 
sound speed, which is often used to define the disc semi-thickness. 

Time dependent evolution of the disc We determine the evolution of the radial 
structure of a non-steady a-disk by solving the diffusion equation (||). To do 
this, we need to use the relation between M s t = 37r(i/)S and S at each radius. 
Interpolation of or piece-wise power law fits to numerical data may be used to 
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Figure 1. Shown is H/r {plot a), £ in units g cm -2 (plot b) and 
T m in K (plot c) vs. r/AU using a logarithmic scale for M s t (in units 
M yr _1 ) = 10" 6 (solid line), 1CT 7 (dotted line), 1(T 8 (short-dashed 
line) and 10~ 9 (long-dashed line) and for a = 1CT 3 . 



represent this relation and more details of these will be published elsewhere. We 
note that they can be used either to compute E from M s t or M s t from E. 

In Figures ||a-b we plot both the curves M s t (S) that we get from the vertical 
structure integrations as described above and those obtained from piece-wise 
power law fits. Figures ||a and ||b are for a = 10~ 2 and 10 -3 , respectively. In 
each case the radius varies between 0.01 and 100 AU. If we calculate M s t using 
the fits with E as an input parameter, the average error is 22, 18, 13% and the 
maximum error is 55, 48, 42% for a = 10- 4 ,10~ 3 , and 10~ 2 , respectively. We 
see that the fits give a good approximation. 

Using the (M s t, E^ relation derived from the integrations, we solved equa- 
tion (^) using explicit finite difference techniques. We considered the situation 
of a disc with initially 0.1 M© for which E oc r _1 , extending to 100 AU. The 
central star had M* = 1 M and for illustrative purposes we adopted a = 10~ 2 . 

In Figure || we show the evolution of E as a function of time. After a time 
~ 10 6 yr, the E profile resembles that of a steady disc in the inner parts justifying 
the assumption of a steady state disc model. After this time, the model is similar 
to that assumed for the primordial solar nebula with E ~ 200 g cm" 2 at 5 AU. 

5. Planetesimal dynamics 

It is thought that planetesimals can be built up from [im sized particles through 
processes of collision, sticking and accumulation occurring in a gaseous medium 
with some degree of turbulence (see the review by Weidenschilling & Cuzzi 
1993). If physical parameters are favourable, particles with a size distribution 
ranging up to ~ a few km can be produced on timescales on the order of 10 yr 
at 1AU. The efficiency of these processes are very uncertain depending as they 
do on sticking probabilities and the degree of particle settling in a turbulent 
medium etc. Planetesimal formation efficiency may be a function of the nebula 
location, being more effective beyond a few AU, where ice has condensed. Here 
we assume that planetesimals with mass m p ~ 10 18 g may form on a sufficiently 
rapid timescale anywhere in the nebula. 

We suppose the planetesimals have number distribution n{m) oc m~ q , for 
some index q. Here the number of planetesimals in the mass range (m, m + dm) 
is n(m)dm. Then most of the mass is distributed in the most/least massive 
objects according as q < 2 or q > 2. The surface density of matter in the form of 
planetesimals is E p . The characteristic number density, N(m p ), for characteristic 
mass m p is then E p /(2m p /i p ). Here h p is the scale height of the distribution 
(~ \/2tt/3 v/Q), where v is the root mean square velocity dispersion. 

For the typical values, E p = 1 g cm" 2 , m p = 10 21 g, and h p = 10 12 cm, the 
mean distance between planetesimals is ~ 10 11 cm. This suggests use of a local 
box model (e.g. Stewart & Wetherill 1988). In this box, the centre of which is in 
circular orbit with the local keplerian angular velocity, planetesimals are assumed 
to move between encounters with constant velocity. In this respect, the effect of 
the central mass is ignored, and the planetesimals are treated using the methods 
of kinetic theory. The idea (Safronov 1966) is that the velocity dispersion of the 
planetesimals is increased by gravitational scattering, enhancing direct collisions 
between them through which they accumulate and grow. 
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Figure 2. M s t in units M Q yr -1 vs. S in units g cm" 2 using a 
logarithmic scale for a = (plot a) and 10~ 3 (plot b). Both the 

curves corresponding to the numerical calculations (solid line) and the 
fits (dashed line) are shown. The label on the curves represents the 
radius, which varies between 0.01 and 100 AU. 
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Figure 3. Solution of the diffusion equation. Shown is £ in units 
g cm" 2 vs. r in AU using a logarithmic scale plotted at times t = 
(solid line), t = 1.36 x 10 4 yr (dotted line), t = 1.1 x 10 5 yr (short- 
dashed line) and t = 1.1 x 10 6 yr (long-dashed line). This run has 
a = 1CT 2 . The total disc mass decreases because of accretion onto the 
central object 
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The local box model fails when the largest planetesimals become isolated. 
That happens when accumulation has reached the stage where there are so few 
that they are in non-overlapping near circular orbits such that they cannot per- 
turb similar mass objects in neighbouring orbits into collision. In this situation 
the central mass can no longer be neglected. To get an approximate idea of 
when isolation occurs at radius r, equate 



This is a statement that there is one object in an annulus of width equal to four 
times its Roche lobe size. Then 



This gives m p = 10 g at 5 AU for E p = 1 g cm~ 2 and M* = 1 M with obvious 
scalings to other surface densities and radii. Thus on the order of 1 may 
be obtained at 5AU. This argument assumes circular orbits which is probably 
reasonable for the largest objects which are circularized through dynamical fric- 
tion. After the isolation stage, planetesimal evolution is probably best followed 
by global N-body methods (e.g. Aarseth et al. 1993). In the gas free case the 
ultimate result is expected to be the formation of terrestrial planets. 

5.1. Gravitational scattering and velocity dispersion 

For a member of a swarm of equal masses interacting gravitationally, the cross 
section for elastic scattering is on average larger than that for a direct collision, 
so long as the velocity dispersion is less than the escape velocity. Thus when 
the velocity dispersion of the swarm is initially small, it will increase because of 
the effect of elastic scatterings. 

The interaction radius for two masses m p to give significant scattering is 
r x = 2Gm p /v 2 . The time between encounters is then t c = l/(Af(m p )7rr 2 Ulog(A)). 
Here the log(A) ~ 10 term accounts for more distant collisions (Binney & 
Tremaine 1987). Using the above, one gets 



For m p = 10 20 g, S p = 1 g cm" 2 and I» = 1 M , one finds at 1 AU that 

t c = 2 x io 17 (V0 4 yr- 



Note that h p is small because the planetesimal dispersion velocities are 
expected to reach the escape velocity for the characteristic mass at a maximum. 
At this point inelastic physical impacts become as important as scattering and 
damp the random motions. Then for m p = 10 20 g, v = 1.7 x 10 3 cm s -1 , 
h p /r = 8 x 10~ 4 and t c ~9x 10 4 yr at 1 AU. 

Thus early planetesimal build up occurs on a rapid timescale. 

5.2. Runaway Accretion 

From the above arguments, the collision time of m p with m p i is inversely propor- 
tional to N(m p ')(m p + m p ') 2 . If this does not decrease with m p i, then collisions 



fn p 



8vr(m p /M*) 1 /3 r 2- 



m,p = 



(8vrS p r 2 ) 3 / 2 (M,)- 1 / 2 . 
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with larger masses are dominant and we expect velocity dispersions to build up 
to the escape velocity of the largest body (Safronov 1966). This occurs for q < 2. 

For q > 2, the largest bodies collide with predominantly smaller ones and 
are circularized by dynamical friction. They can then accrete efficiently from 
the smaller ones which move with a velocity dispersion small compared to the 
escape velocity from them. 

The accretion rate from planetesimals with mass m p i and velocity dispersion 
v p i is enhanced by gravitational focusing. Thus the growth of the largest mass 
m p say, with radius R p is given by 

dm v AT . . _, 2 _ / 2Gm p \ 
— = N(m p ,)m p ,nR p v p , ^1 + _ j . 

Runaway is caused both by the increase of R p with m p and the gravitational 
focusing term in brackets. The timescale for growth to isolation mass is compa- 
rable to the encounter time, t c , indicated above. But note that growth may be 
slowed down through the effect of encounters between neighbouring runaways 
producing an increase in the velocity dispersion that propagates to all compo- 
nents of the system (Ida & Makino 1993) 



6. Giant planet formation 

After a solid protoplanetary core grows to a critical mass of around a few M®, 
the surrounding gaseous atmosphere can no longer grow quasi-statically in mass 
along with it. A process of collapse ensues possibly leading to dynamical ac- 
cretion and mass growth to values characteristic of giant planets. We review 
the theory of this below. Such a critical core mass model is supported by the 
indication from models of Jupiter that it has a solid core of ~ 5-15 M® (Podolak 
et al. 1993). 

6.1. Basic equations governing a protoplanetary envelope 

Let w be the spherical polar radius in a frame with origin at the centre of 
the planet's core. We neglect the rotation of the planet around both its own 
spin axis and the disc spin axis. We assume that the envelope is in hydrostatic 
equilibrium and spherically symmetric, so that: 

^ = -gp, (18) 

aw 

where g = GM/w 2 is the acceleration due to gravity, M(w) being the mass 
contained in the sphere of radius w (this includes the core mass if w is larger 
than the core radius). Mass conservation gives: 

dM 9 

— = 4irm 2 p. (19) 
dzu 

The thermodynamic variables in the protoplanet envelope are such that the 
equation of state of an ideal gas does not normally apply. Here we adopt the 
state-of-the-art Chabrier et al. (1992) equation of state for a hydrogen and 
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helium mixture. We fix the abundances of hydrogen and helium to be 0.7 and 
0.28 respectively. 

The equation of radiative transport is: 

dT_ -3np L 

dm 16ar 3 Attw 2 ' 1 ' 

where L is the luminosity carried by radiation. Denoting the radiative and 
adiabatic temperature gradients by V ra d and V a d respectively, we have 

_/91nT\ 3nL core P 
Wrad ~ \d\nP) rad ~ MvaGMT^ 

and 

Mis?)/ < 22 > 

with the subscript s meaning that the derivative has to be evaluated at constant 
entropy. 

We assume that the only energy source comes from the core which outputs 
the core luminosity L core , given by: 

j - m GM core 

'-'core — lvl core i y^ ) 

r core 

where M core and r core are respectively the mass and the radius of the core, 
and Mcore is the rate of accretion of planetesimals onto the core. We note 
that it is customary to take, instead of L core , the luminosity supplied by the 
gravitational energy which the planetesimals entering the planet atmosphere 
release near the surface of the core (see, e.g., Mizuno 1980; Bodenheimer &: 
Pollack 1986). However, when the mass of the atmosphere is small compared 
to that of the core, these two luminosities are comparable, providing we take 
for M core the rate of accretion of planetesimals onto the core during the phase 
of accretion of the atmosphere and not during the phase of the core formation. 
These two rates may be different if the core has migrated in the disk before 
accreting the atmosphere. 

If V ra d < V 'ad, t ne medium is convectively stable and the energy is trans- 
ported only by radiation. In that case L = L core . 

When V ra d > V a d, some energy is transported by convection. In that 
case, L core = L + Lconv, where L con v is the convective luminosity. We use the 
expression for L CO nv given by mixing length theory (Cox & Giuli 1968): 



dT \ _ f dT 
dzu J , \ dm 




dT) 



p 



(24) 



where A m i = \a m iP/(dP/dm)\ is the mixing length, a m i being a constant of 
order unity, (dT/dm) s = V a dT (din P/ dm), and the subscript P means that the 
derivative has to be evaluated for a constant pressure. The quantities (dp/dT) p 
and V a d are given by Chabrier et al. (1992). 
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6.2. Boundary conditions 

We suppose that the planet core has a uniform mass density p core , here taken to 
be 3.2 g cm -3 . The core radius, which is the inner boundary of the atmosphere, 
is then given by: 

' . (25) 



v 47T/3 core , 

The outer boundary of the atmosphere is taken to be at the Roche lobe 
radius of the planet: 

2 f M pi\ 1/3 (0fi , 
rL = 3{3Mj r ' (26) 

where M p i = M core + M atm is the planet mass, M atm being the mass of the 
atmosphere, and r is the location of the planet in the disk (i.e. the separation 
between the planet and the central star). 

To avoid confusion, we will denote the disk mid-plane temperature, pressure 
and mass density at the distance r from the central star by T mic i, P m id and Pmid-, 
respectively. 

At w = tl, the mass is equal to M p i, the pressure is equal to P m id and the 
temperature is given by: 



T 4 T L^core i 

mid+ A-narl) ' [ ' 



1/4 



where 



T L — {Pmidi Tmid) Pmid r L- 



The condition at Til — fcore is that the mass is equal to J\d core there. 
6.3. Model calculations 

At a given disk radius r and for a given core mass J\d core , we solve the equa- 



tions (jig), ( |19|) and fl20D with the boundary conditions described above to get 
the structure of the envelope. The opacity law adopted was the same as that for 
the disk models. We note that when the density gets large, the interior of the 
envelope becomes convective so that the value of the opacity does not matter 
there. 

The equations are integrated using the fifth-order Runge-Kutta method with 
adaptive step-size control (Press et al. 1992). We guess a starting value of M atm 
and integrate the equations from w = ri down to the core surface w = r core . 
We then iterate the integration, adjusting M atm at each step, until the solution 
gives M = M core at w = r core with some accuracy. 

At each radius r, for a fixed M core , there is a critical core mass M cr n (which 
increases as M core increases) above which no solution can be found, i.e. there 
can be no atmosphere in hydrostatic and thermal equilibrium confined between 
the radii r core and tl around cores with mass larger than M cr , t . This is because 
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when the core mass is too large, the atmosphere has to collapse onto the core in 
order to supply adequate luminosity to support itself. For masses below M cr i t , 
there are (at least) two solutions, corresponding to a low-mass and a high-mass 
envelope respectively. 

In Figure ^ we plot M v \ versus M core for different M core (between 10~ 6 
and 10 -11 M yr _1 ) at a radius of 5 AU and for T mic i = 140.05 K and P m id = 
0.13 dyn cm~ 2 . These values of the temperature and pressure are obtained 
from the vertical structure integrations described above when the parameters 
a = 10~ 2 and M st = 10~ 7 M yr -1 are used at r = 5 AU. 

The critical core mass, which decreases when M core _ decreases, is found to 
be 22.5 M e for M core = 10" 6 M e yr" 1 and 2.5 M e for M core = 10 -11 M e yr" 1 . 
These values are slightly larger than those found by Bodenheimer & Pollack 
(1986). The difference may be accounted for by the fact that we do not calculate 
the luminosity in exactly the same way as they do. Also we use a slightly different 
boundary condition for the temperature at the surface of the planet. 

7. Disc protoplanet interaction 

Once the planetary mass has attained values of around an earth mass or higher, 
dynamical interactions with the surrounding disc matter become important, 
leading to phenomena such as inward orbital migration and gap formation (Lin 
& Papaloizou 1993; Ward 1997; Lin et al. 1998). 

Korycansky & Papaloizou (1996) considered the perturbed disc flow around 
an embedded protoplanet when the imposed viscosity v = 0. They used a 
shearing sheet approximation in which a patch, centered on the planet, coro- 
tating with its orbit is considered in a 2D approximation. For unit of length 
r t = r(MpifM^) 1 / 3 was adopted, where r is the planet's orbital radius. We re- 
mark that rt = (3 4 / 3 r^)/2, is a multiple of the Roche lobe radius used above. 
When the basic equations are expressed in dimensionless units, the only param- 
eter defining the problem is Ai = rt/(c s /Q), being essentially the ratio of Roche 
lobe radius to disc semi-thickness. 

The velocity v viewed in the rotating frame was split into components 
involving a vortical function T and a potential <3?: 

v = Vxir + V$ (28) 

For a steady flow, there is also a stream function ^f, such that 

v = (V x z^>) /£ (29) 

Parameters relating to the flow for A4 = 0.7 are plotted in figures 5 and 6. Here 
Cartesian coordinates are adopted with origin at the centre of the protoplanet 
and x axis pointing along the line joining the central star to the protoplanet. 
The perturbed surface density (figure 5) shows trailing shock waves behind which 
there is a strong surface density enhancement giving rise to pronounced wakes. 
Also notable is the prograde disc flow around the protoplanet (figure 6). Dissi- 
pation in the shock waves, which become strong once A4 > 1, eventually results 
in gap formation (Lin & Papaloizou 1993). 
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Figure 4. M v \ in M ffi versus M core in M ffi for M core between 10 6 
and l(T n M e yr" 1 at a radius r = 5 AU and for T mid = 140.05 K and 
Pmid = 0.13 dyn cm -2 . From one curve to another, starting from the 
right, M core decreases by a factor 10. The critical core mass increases 
with Mcore-, varying between ~ 2 and 22.5 M®. 
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Figure 6. Velocity field (including background) near the protoplanet 
taken from Korycansky & Papaloizou (1996). Left: Flow field. Right: 
Velocity v y , along the line y = 0. 
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Protoplanet-disc interaction leads to gap formation and orbital migration, 
which, together with tidal interaction with the central star (Terquem et al. 
1998) can lead to massive planets in circular orbits, with periods of a few days, 
as observed (Butler et al. 1997). The reader is referred to the article by Lin et 
al. (1998) for an account of the dynamical phenomena which can play a role in 
determining the final orbital configuration of planetary systems. 
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